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Structural  mechanics,  finite  elements,  nonl inear  equations , 
iterative  solution  methods 


'Tire  solution  of  nonlinear,  transient  finite  element  prob¬ 
lems  may  be  achieved  by  using  step-by-step  integration  of  the 
equations  of  motion  combined  with  a  Newton  solution  of  the 
resulting  nonlinear  algebraic  equations.  The  use  ot  Newton 
type  methods  leads  to  a  set  of  linear  simultaneous  algebraic 
equations  whoso  solution  gives  the  next  iterate,  for  very. 
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large  problems  the  solution  of  the  large  set  of  linearized 
equations  may  be  a  formidable  task  -  often  consuming  more  than 
half  of  the  computing  effort  when  performed  by  a  direct  method 
based  upon  Gauss  elimination.  Accordingly,  it  is  of  consider¬ 
able  importance  to  investigate  alternative  methods  to  solve  the 
problem.  The  present  study  presents  results  obtained  by  using 
•a  Preconditioned  Conjugate  Gradient  Method  (PCG)  described  in 
[7]  and  a  Preconditioned  Lanczos  Method  (PLM)  described  in  [6] 
to  solve  a  variety  of  numerical  examples.  Based  upon  results 
obtained  it  is  evident  that  a  significant  reduction  in  overall 
effort,  compared  to  direct  solutions,  may  be  achieved  using  the 
preconditioned  methods. 
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t.  Introduction 


The  finite  element  method  of  discretization  is  used  to  reduce  many  complex  continuum 


problems  to  discrete  systems.  Although  this  reduction  is  the  most  important  step  in  the  overall 


analysis  ol  a  structure,  solving  the  resulting  discrete  problem  is  often  far  from  trivial.  In  gen¬ 


eral,  the  reduced  system  is  nonlinear  and  an  iterative  method  must  be  employed  to  arrive  at  the 


solution.  Most  solution  methods  are  based  on  some  form  of  Newton's  method  in  which  the 


nonlinear  problem  is  linearized,  using  an  initial  approximation,  to  arrive  at  a  linear  set  of  simul¬ 


taneous  algebraic  equations.  The  solution  of  the  set  of  linear  equations  leads  to  a  correction  ot 


the  initial  approximation.  When  solving  the  linear  equations,  one  should  not  loose  sight  of  the 


primary  objective,  solving  the  nonlinear  problem. 


Iterative  methods,  such  as  the  conjugate  gradient  or  Lanczos  method,  are  among  the 


many  methods  that  may  be  used  to  solve  systems  of  linear  equations.  The  advantage  of  these 


methods,  when  used  as  the  inner  loop  of  the  Newton  iteration,  is  twofold. 


(i)  i  he  linear  equation  may  be  solved  to  any  desired  level  of  accuracy  as  governed  by  the 


Newton  iteration. 


(ii)  A  considerable  reduction  in  storage  can  be  achieved  when  no  triangular  factorization  need 


be  performed. 


In  1 6]  a  method  was  developed,  based  on  the  preconditioned  Lanczos  method,  to  realize 


some  of  the  advantages  of  iterative  methods.  In  this  previous  study,  the  triangular  factors  of 


the  initial  tangent  matrix  were  used  to  form  a  preconditioning  matrix  for  the  subsequent  solu¬ 


tion  steps  In  the  present  study  we  have  eliminated  factorizations  by  employing  other  precondi¬ 


tioners  and  further,  have  reduced  the  storage  needs  of  the  method  Accession  For  j 
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2.  A  Preconditioned  Conjugate  (iradienl  Method 

An  essential  step  in  nonlinear  analysis  of  structures  using  Newton’s  method  (or  a  variant 
such  as  modified  Newton  or  quasi-Newton  methods)  is  solving  a  linear  system  ol  algebraic 
equations.  The  preconditioned  conjugate  gradient  method  (hereafter  called  PC’(i)  is  one  of  the 
many  procedures  for  solving 

r  =  b  Ax  =  0  t2  1  < 

where  A  is  an  n>  n  symmetric  positive  definite  matrix  (which  for  finite  element  calculations  is 

sparsely  populated)  and  b  is  the  right-hand  side  vector.  In  the  case  of  static  analysis.  A  is  the 

current  tangent  matrix  and  in  the  case  of  dynamic  analysis,  A  depends  on  the  mass,  damping 

and  tangent  stiffness  matrices,  as  well  as  the  time  increment. 

The  initial  popularity  of  the  conjugate  gradient  method  was  due  to  a  number  of  factors. 
In  exact  arithmetic  the  method  required  a  maximum  of  n  iterations  to  solve  (2.1)  which  made 
the  method  superior  to  other  iterative  methods.  In  fact  conjugate  gradient  is  in  the  class  of 
semi-iterative  methods  which  also  includes  the  Lanczos  algorithm  llOl  The  disadvantage  ol 
direct  methods  is  their  large  storage  demands  for  keeping  the  factors  of  A.  The  only  interface 
between  the  conjugate  gradient  method  and  A  is  through  the  product  Av  for  a  given  vector  v. 
This  is  an  elegant  way  of  taking  advantage  of  sparsity  of  A  which  has  the  added  advantage  that 
A  need  not  be  known  explicitly  but  only  a  means  of  computing  the  matrix-vector  product  is 
required 

The  popularity  of  the  conjugate  gradient  method  vanished  once  it  was  found  that  under 
certain  conditions  the  method  required  as  many  as  5«  or  bn  steps  to  reduce  the  residual  to  the 
desired  level.  I  his  degradation  is  due  to  the  strong  influence  of  round-off  error. 

The  addition  ot  preconditioning  eliminated  this  difficulty  Instead  of  solving  (2  I)  we 

solve 


P  Ax  -  P  'b  (2  2) 

tor  some  appropriate  choice  of  P  The  object  then  is  to  choose  P  such  that  the  coefficient 
matrix  of  (2  2>  is  well  conditioned 


sfo  *.*.  e'</ 
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Theoretical  considerations  suggest  that  at  the  end  of  each  iteration  of  C (i  the  residual 
^  ■  1 

norm  is  reduced  bv  a  factor  — -  when  solving  (2.1)  where  x  <s  the  condition  number  of  A. 

•Jk  +  I 

defined  by  k  =  it  A  i!  ||  A  !||.  See  ll]  for  more  details.  Note  that  when  k  =  1,  one  iteration  is 
sufficient  to  solve  the  equation.  This  provides  us  with  a  guideline  for  choosing  P  l  or  a  well 
chosen  P  only  a  few  iterations  reduce  the  residual  norm  to  the  desired  level  Mere  we  gi\e  an 
outline  for  the  preconditioned  conjugate  algorithm: 

Given  an  initial  guess  x„,  a  positive  definite  preconditioning  matrix  P.  the  matrix  A  and  the 
right  hand  side  b: 

( 1 )  Set  p,,  r„  b  Ax,) 

(2)  Solve  Pd, i  -  r, .  for  d,, 

(4)  for  A  -  0.  I.  2,  until  convergence  do 

(at  <>.  1  <r,  ,d„  )/(pA,Ap* ) 

>b)  X;.(  -  x<  4 
<c>  r i .  |  =•  r.  -(t(  Ap, 

(d)  Solve  I’d t .  |  -  rA ,  | 

(e)  -  <r.  .  ,.d* .  | )/(rA ,d* ) 

<f)  p.i  =  dk.,  \iik p* 

The  operation  (v,u)  denotes  the  inner  product  v'u.  This  algorithm  generates  a  sequence 
of  approximations  to  the  solution  x  with  a  corresponding  residual  vector  rA.  The  termination 
criterion  can  be  chosen  based  on  these  quantities.  In  addition  to  storage  demands  for  A  and  P 
the  algorithm  requires  storage  for  4  vectors.  The  total  number  of  operation  per  iteration  is 
V/l  t  S/P  t  5 A.  where  X/.1  and  M/P  are  the  number  of  operations  for  forming  An  and 
P  1  v  for  a  given  ii  and  \ 
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3.  Splitting  Methods 

Next  we  turn  *o  a  topic  which  at  first  sight  may  seem  unrelated  to  the  solution  of  non¬ 
linear  algebraic  equations.  Consider  the  system  of  first  order  differential  equations 

x  =  f(x,/)  (3. 1 ) 

where  x  is  an  n-dimensional  vector,  the  superposed  dot,  (  ),  denotes  differentiation  with 

respect  to  time  and  f  is  a  function  of  the  unknown  vector  x  and  t. 

We  consider  a  special  form  of  f  which  can  be  written  as  a  sum  of  its  subcomponents  f 

f-  ff,  (3.2) 

,  i 

Under  these  conditions  the  original  problem  can  be  thought  as  a  sum  of  s  subproblems 

x  =  f,(x,r)  /  =  1 _ _  .v  (3.3) 

In  the  case  of  finite  element  discretization  of  the  spatial  domain  the  sum  in  (3.2)  ranges  over 

the  elements  or  a  set  of  elements.  In  other  cases  the  splitting  may  be  formed  by  other  means. 

one  of  which  is  demonstrated  in  the  following  section. 

consistent  algorithm  for  the  solution  of  (3.1),  based  on  the  notion  of  a  splitting  tech¬ 
nique  {21,  can  now  be  constructed  as  a  product  of  algorithms  for  the  sub-problems.  In  other 
words,  write  the  algorithm  for  (3.3)  as 

0  4) 

where  .S'1'"  is  an  operator  denoting  the  algorithm  and  the  index  m  ranges  over  the  increment  in 
time.  h.  Then  the  algorithm  for  (3.1)  can  be  written  as 

x,,,.]  =  .V’Mx,,,]  (3.5) 

where 

.V"1  =  (I  .V,'"’  (3.M 

,  i 

One  of  the  disadvantages  of  the  splitting  method  is  its  low  accuracy.  The  best  that  these 

methods  can  achieve  is  second  order  accuracy.  That  is  the  truncation  error  is  of  the  order  of  hx 

at  best  In  the  sequel  wc  will  use  the  above  procedure  to  construct  a  preconditioning  matrix  for 
the  conjugate  gradient  algorithm  described  in  section  2.  The  inherent  inaccuracy  of  the  splitting 


k  • 

l- 


r  ■' 

k  . 


method  poses  no  problem  since  the  algorithm  is  used  only  as  a  preconditioner  and  therefore 
one  can  obtain  very  high  accuracies  through  the  conjugate  gradient  iteration. 


4.  Solution  of  Static  Problems 

Consider  ihe  system  of  linear  first  order  differential  equations 

rx  +  Ax  =  b  (4  1' 

where  r  is  a  given  parameter.  Formally  the  solution  to  equation  (4.1)  is 

x(t)  -  <■  Al  Mx„  -  A  !b)  +  A  'b  (4.2) 

where  x„  =  x(0),  is  an  initial  condition.  We  observe  from  (4.2)  thet  as  t  ends  to  infinity  xl/> 

converges  to  the  solution  of  (2.1)  for  t  >  0.  Consiquently  (4.1)  may  be  utilized  to  solve  the 
linear  equations  (2.1).  Indeed  this  approach  has  been  suggested  previously  (e  g.,  see  [9])  in 
general  the  exponential  of  a  large  matrix  cannot  be  easily  computed  and  a  numerical  solu>  )f 
(4  1)  must  be  used  In  order  to  achieve  a  soluion  of  (2.1)  the  numerical  solution  to  U  must 
be  assymptotically  correct  for  infinite  h ,  or  a  very  large  number  of  time  steps  must  1  us  '  to 
compute  the  solution  at  infinite  time.  Here  we  are  not  concerned  with  constructing  an  accurate 
solution  to  (4.1),  rather  we  consider  the  method  as  a  means  of  constructing  a  suitable  precondi¬ 
tioning  matrix  for  the  conjugate  gradient  algorithm  described  above. 

Splitting  methods  may  be  applied  to  any  problem  of  the  form 

x  »  Bx  (4.3) 

where  B  is  an  additive  operator  defined  by 

B  =  £  B,  (4.4) 

,  -  i 

such  that  the  equations 

x  =  B,x  i  =  l,...,s  (4.5) 

are  significantly  easier  to  solve  than  the  original  equations.  The  time  stepping  algorithm  for  the 

global  problem  is  then  the  product  of  all  the  time  stepping  algorithms  for  the  subproblems  with 
a  fractional  time  step  h/s  (2). 

The  coefficient  matrix  A  in  (2.1)  may  be  written  as  the  sum  of  its  diagonal  matrix.  I),  a 
strictly  lower  triangular  matrix,  L,  and  a  srictly  upper  triangular  matrix  such  that 

A  =  C/?D  +  L)  +  C/?D  +  I„) 1  (4  6) 

The  associated  subprobiems,  x  =  -O/2D  +■  L)x  and  x  =  -C/.D  ■*  I. ;  )x  can  be  solved  easily 
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Applying  a  backward  difference  method  with  a  time  step  h! 2  to  each  of  the  subproblems,  we 
arrive  at 


I  4  /(D  +  21.) 

1 

I  4  -^-(D  4  2L 7  ) 

1 

h  . 

---  b  t-  \ 

4r 

4t 

It 

where  x,„  is  an  approximation  to  xirtih).  For  an  initial  condition  x0  =  0  we  get  an  approxima¬ 


tion  to  x(  /;) 


[  Jf 

1  4  --  (D  4  2L) 

1 

I  4  /'  (D  4  2L7  ) 

1 

I2t 

4r 

4t 

which  is  compared  to  the  exact  solution 


x(h)  -  (l  -  e  A'  ]a  1  b  (4.9) 

Comparing  equations  (4.8)  and  (4.9)  suggests  that  the  coefficient  matrix  in  (4.8)  may  be  a  good 

approximation  to  A  1  for  large  h,  and  may  therefore  be  an  effective  preconditioning  matrix. 

The  scalar  factor  may  be  ignored  for  preconditioning  purposes. 

It 

When  using  (4.8)  in  conjunction  with  the  conjugate  gradient  algorithm  of  section  2  the 
preconditioning  matrix  becomes 

P  =  (I  +  <o/2D  4  wLHI  4  <o/ 2D  +  toL')  (4.10) 

where  tu  -  /;/2r  is  now  a  free  parameter. 

To  simplify  the  choice  of  w  we  scale  the  stiffness  matrix  A  such  that  diagonal  of  A  is 
unity.  The  resulting  matrix  is  A  =  D  AD  The  system  of  equations  (2.1)  now  becomes 

A  x  -  b  (4  11) 

where  x  D  x  and  b  I)  b 

I  he  preconditioned  matrix  must  now  be  applied  to  (4  11)  resulting  in 

Pdt  toL) (I  4  „[')  U  12i 

where  A  I.  f  I.  It  is  easy  to  show  that  preconditioning  (4.11)  using  P  is  equivalent  to 

preconditioning  (2  1)  with 


P  (D  +  tolJD  1  (D  -f  toL') 


(4.13) 


This  can  be  identified  as  a  member  of  the  class  of  incomplete  Choleski  preconditioners  |3| 


Note  that  when  «j  =  0,  P  becomes  the  diagonal  matrix  D,  resulting  in  the  simplest  form  o I' 
preconditioning;  diagonal  scaling.  When  <»  I  then  P  =  A  f  LI)  'I.'  where  we  note  that  the 
error  matrix  LD  'L'  is  rank  deficient  since  L  has  zero  diagonals.  If  the  norm  of  D  is  largei 
then  the  norm  of  L  then  the  norm  of  the  error  matrix  will  be  small  compared  to  the  norm  of 
A.  consequently,  for  most  problems  it  is  expected  that  the  optimum  w  will  be  close  to  unity 


5.  Solution  of  Dynamic  Problems 

We  next  construct  u  preconditioning  matrix  for  the  linear  system  of  equations  arising  in  a 
step-hv-step  algorithm  tor  dynamic  analysis  of  linear  and  nonlinear  structures.  In  particular,  sc 
consider  the  Newmark  algorithm  and  the  preconditioning  matrix  follows  from  the  splitting 
method  of  section  3,  in  much  the  same  way  as  for  the  static  problem. 

Consider  the  linearized  equations  of  motion 

Mii  f  ku  =  f  <  s  1 ) 

where  M  is  the  diagonal  mass  matrix,  K  is  the  stiffncs.,  matrix,  f  is  the  external  load  vector  and 

u  is  the  response  ol  the  structure  lor  simplicity,  we  ignore  damping  effects;  however,  all  of 

the  following  results  may  he  extended  easily  to  the  damped  case.  Accordingly,  the  linearized 

system  of  equations  arising  at  every  time  step  of  the  Newmark  method  is 


where 


A  x  -  b 


(5.2) 


and 


K  4 


■ TM 

0A  i2 


(5  3) 


b  =•-  f,  .  a/  4  -  '  Miu,  4  A/v,  4  C/j  -  /HArY  J  (5  4) 

/fAr 

Here  v  and  a  are  velocity  and  acceleration  vectors,  respectively,  A r  is  the  specified  time  incre¬ 
ment,  r  is  the  time  and  x  is  now  the  increment  of  displacement  response.  The  Newmark 
parameters  are  chosen  such  that  /f  ^  {'/:  4  y):/4  with  y  Si  '/?  which  ensures  unconditional 
local  stability.  The  discretizations  in  lime  are 

u,  -  u.  *  Arv.4  i-'k  Ar  [(1  -  2/J)a,  4  2/3a,  ,  a,)]  (5.5) 

C  a,  v  +  A / ( I  y)a,  4  yAra,  .  a,  t5.M 

I  he  object  is  to  solve  (5.2)  without  forming  the  factors  of  A. 

A  splitting  method  similar  to  the  one  used  for  equation  <4.1 )  can  now  be  applied  to  equa¬ 
tion  (51)  The  matrix  resulting  front  the  splitting  algorithm  can  then  be  used  as  a  precondi¬ 


tioner  for  1 5.2)  (  onsidei 


<5  7) 


P  =  (I.  +  -  -  IM)M  '(1/  4  —  - M) 

pAr  pAt2 

where  K  ---•  L  •  I. Multiplying  out  the  terms  in  (5.7),  we  obtain 

P  -  --  (/iArLM  'L'  +  L  +  L'  +  M) 

pAr  pAt2 

=  \pArlM  1 L '  4  A I 

pAr 

-  lE(Ar)  +  Al  (5.8i 

PAr 

where  E<A r)  -  pArLM  ‘L 7 . 

The  preconditioned  conjugate  gradient  algorithm  of  section  2  is  invariant  under  the  scaling 
of  the  preconditioning  matrix;  therefore,  (5.8)  shows  that  P  will  tend  quadratically  to  the 
dynamic  stiffness  matrix  A  as  the  time  step  diminishes.  In  other  words,  E  tends  to  the  zero 
matrix  quadratically  in  At  We  see  later  that  this  characteristic  results  in  an  effective  precondi¬ 
tioning  and  the  solution  of  equation  (5.2)  is  obtained  in  very  few  iterations  of  the  precondi¬ 


tioned  conjugate  gradient  algorithm. 


6.  l.anczos  Algorithm  for  Solution  of  the  Linearized  Problem 

The  discretization  of  nonlinear  structural  mechanics  problems  and  linearization  ol  the 
resulting  nonlinear  algebraic  equations  for  application  of  a  Newton  type  methods  normal!;,  will 
lead  to  a  symmetric  system  of  linear  algebraic  equations,  (2.1). 

In  Section  2  we  described  the  use  of  the  conjugate  gradient  method  to  solve  the  linear 
system  of  equations.  In  this  section  an  alternative  solution  technique,  a  simple  l.anczos 
method,  is  presented 

Iterative  methods  often  have  been  used  in  numerical  analysis  for  the  solution  of  large  sys¬ 
tems  of  equations  1  he  C  onjugate  Gradient  method  is  one  such  technique  introduced  in  l^S2 
by  Mestenes  and  Stiefel  14).  In  the  same  year  Lanczos  published  his  method  of  minimized 
iteration  which  wds  initially  introduced  lor  computing  the  eigen  pairs  of  a  large  symmetric 
matrix.  1  aiic/os  and  Householder  |5|  pointed  out  the  intimate  connection  between  the  two 
approaches  These  methods  have  several  attractive  features  in  common  There  is  no  need  for 
A  to  have  further  special  properties,  such  as  banded  form,  no  acceleration  parameters  have  to 
be  estimated,  and  the  storage  requirements  are  only  a  few  n-vectors  in  addition  to  the  storage 
needs  ol  A 

6.1.  The  I  anezos  Algorithm 

1  or  certain  applications  of  the  finite  element  method,  especially  in  nonlinear  problems,  u 
is  usual  to  have  on  hand  an  initial  approximation  x1'  to  the  true  solution  of  2.1.  The  problem  is 
now  to  find  a  correction  x  to  be  added  to  x".  Then 

Ax  =  r„  (6  I ) 

where 

r„  -  b  Ax1'  U\2> 

The  Lanc/.os  algorithm  may  be  described  very  simply  as  a  process  of  constructing  the  weak 
lorm  ol  equation  2.1  from  a  very  special  subspacc.  The  subspace  under  consideration  is  gen¬ 
erated  from  the  set  of  i  vectors  (  r„  ,  Ar„  ,  •  ,  A'  1  r0  known  as  the  Krylov  suhsnaco  |h>| 
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To  construct  the  weak  form  it  would  be  simpler  if  an  orthonormal  set  of  vectors,  say 
(q,,q:,  ,q,>,  "ere  available.  This  can  be  achieved  by  applying  Gram-Schmidt  orthogonaliza- 

tion  to  the  Krylov  vectors.  Initially,  this  appears  to  be  an  expensive  way  of  obtaining  an  ortho¬ 
normal  base  vectors,  however  this  process  can  be  simplified  when  the  following  two  facts  are 
used  [10]: 

<i>  The  use  of  Aq,  and  A  Tn,  for  orthogonalization  against  the  previous  q  vectors  and  normal¬ 
ization  of  the  resulting  vector,  leads  to  the  same  vector  q,. 

(n)  The  vector  Aq,  is  orthogonal  to  q(,q>,  ,q,  j. 

Consequently  it  is  sufficient  to  orthogonalize  Aq,  against  q,  |  and  q,  to  obtain  the  next  orthog¬ 
onal  vector  Accordingly. 

r,  =  d. .  i q , ■  i  =  Aq,  -  or,q/  -  (i , q,  ,  (h.3) 

where  o  =  q  /Aq  and  ji  q/  |Aq(.  It  is  important  to  note  that  the  vectors  (qi ,q^,  •  .q  >) 
are  not  needed  in  equation  6.3  .  This  defines  one  step  of  the  simple  Lanczos  algorithm  The 
normalization  of  r,  results  in  q, ,  It  is  easy  to  show,  by  looking  at  q/,  jr,,  that  fi,< ,  =  ||r.  !|. 

The  special  choice  for  the  base  vectors  of  the  subspace  has  an  additional  advantage.  The 
projection  of  A  onto  this  subspace  is  a  tridiagonal  matrix,  T,. 

«,  Hi 
(h  "2  Hi 

T ,  =  Q /  AQ,  =  ■  (6.4i 

«/  i  H  i 

H,  <*, 

where  the  q  vectors  form  the  columns  of  the  matrix  Q,,  Q,  (q|,q:,  •  •  ,q,L  This  fact  was 
realized  soon  after  l  anczos  introduced  his  method  and  the  algorithm  was  put  to  use  as  a  pro¬ 
cess  tor  the  orthogonal  transformation  of  a  matrix  to  tridiagonal  form.  Despite  its  additional 
attractions,  the  Lanc/os  process  gave  way  to  Givens'  method  in  1954  and  lain  to 
Householder's  method  in  1958 

The  relationships  that  define  the  simple  Lanczos  algorithm  can  now  be  summari/eii  in  the 


following  three  equations. 
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Q/Q,  -  1, 

AQ,  -  Q,T (  =  r  e/ 

Q/r  =  0 


(6  5.al 
* 6  5.br 
16. 5.c) 


where  e,  is  the  /-th  column  of  the  jx j  identity  matrix  l(.  Setting  q(l  =  0  and  using  r  as  the 
starting  vector,  the  Lanczos  algorithm  may  then  be  described  as 
Given  r,,.  set  yd,  -  ||r0||,  for  /  =  1,2,  •  •  repeat 

«"  -■ -j,1 

(2)  u,  —  Aq, 

(3)  r.  —  u  -  /i,q,  | 

(4)  a,  —  q/r, 

(5)  r,  <—  r  -  «,q, 

<b)  Hr\  —  Hr, II 

While  a  direct  use  of  the  simple  Lanczos  algorithm  usually  leads  to  numerical  difficulties, 
thus  requiring  some  reorthogonalization  of  the  vectors,  with  care  in  selecting  the  precondition¬ 
ing  matrix  these  difficulties  are  avoided.  In  our  test  of  the  algorithm  to  date  no  difficulties  have 
been  encountered  in  using  the  simple  preconditioned  Lanczos  algorithm  (see  Numerical  Lxam- 
ples). 
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1 .  Modification  to  the  FEAP  Program. 

The  finite  element  computer  program  FEAP  is  described  in  Chapter  24  of  Zienkiewicz 
(13)  and  forms  the  basis  for  all  computations  performed  as  part  of  the  current  effort.  Several 
basic  modifications  have  occurred  since  the  original  program  was  published.  These  include:  (a) 
replacement  of  the  original  subprograms  ACTCOL  and  UACTCL  for  the  direct  solution  of 
linear  algebraic  equations  for  symmetric  and  unsymmetric  problems,  respectively,  by  the  single 
solver  DASOL  which  is  probably  the  most  efficient  implementation  of  a  variable  band,  active 
column  equation  solver  available  today  for  virtual  memory  machines;  (b)  inclusion  of  the  New- 
mark  method  to  integrate  the  equations  of  motion  for  linear  and  nonlinear  problems  in  struc¬ 
tural  mechanics,  (c >  the  ability  to  construct  a  tangent  stiffness  matrix  and  a  residual  force 
simultaineously  (instead  of  using  TANG  followed  by  FORM,  e.g.,  see  below);  (d)  a  conver¬ 
gence  criteria  based  upon  the  current  increment  in  energy  (e.g..  Ax7  r);  and  (e)  inclusion  of  a 
general  data  storage  structure  for  nonlinear  materials  which  require  additional  information  to 
the  displacements  in  order  to  evaluate  the  current  stress  state  at  each  integration  point  in  an 
element 

In  the  work  reported  here  the  program  FEAP  has  been  extended  further  to  include  the 
capabi.  of  solving  problems  using  the  Preconditioned  Conjugate  Gradient  Method  (PCG)  or 
the  Preconditioned  Lanczos  Method  (PLM)  described  above.  In  addition,  a  line  search  algo¬ 
rithm  has  been  incorporated  for  use  with  any  of  the  solution  methods  -  i.e.,  direct,  PCG,  or 
PL.M  As  noted  above,  the  preconditioning  matrices  require  a  knowledge  of  the  nonzero  terms 
in  the  global  tangent  stiffness  matrix  (as  well  as  the  mass  matrix  for  dynamics  problems).  Since 
it  is  not  necessary  to  factor  the  preconditioning  matrix,  which  would  cause  fill  in  the  nonzero 
structure,  we  have  developed  a  direct  means  to  construct  the  array  containing  the  nonzero 
terms. 

The  algorithm  to  construct  the  locations  of  the  nonzero  terms  in  the  compressed  tangent 
stiffness  matrix  may  be  summarized  by  the  following  steps: 
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1  )  Make  a  list  of  the  elements  which  are  attached  to  each  node.  In  FEAP  this  step  is  accom¬ 
plished  by  constructing  dynamically  dimensioned  arrays  which  will  contain  only  the  terms 
associated  with  the  actual  nonzero  structure;  so  that  no  storage  will  be  wasted.  Accord¬ 
ingly,  the  first  step  is  divided  into  two  parts  -  the  first  to  determine  dimensioning  and  the 
second  to  obtain  the  actual  elements  connected  to  each  node. 

2.)  For  each  nodal  degree-of-freedom  use  the  list  of  elements  to  find  all  the  other  nodal 
degree-of-freedoms  which  are  connected.  Since  we  are  currently  considering  symmetric 
equations  only  the  terms  above  the  diagonal  entry  are  constructed. 

The  above  two  steps  have  been  incorporated  into  the  set  of  subprograms  which  are  listed 
in  Appendix  A.  In  FEAP  the  array  of  nonzero  terms  may  be  constructed  by  using  the  macro 
command  CTAN  (compressed  tangent  stiffness  matrix)  instead  of  the  usual  TANG  (note  that 
there  is  no  equivalent  compressed  array  for  UTAN  since  neither  of  the  preconditioned  solvers 
is  programmed  to  handle  unsymmetric  equations).  Thus  a  typical  Newton  iteration  using  the 
PCG  a'gorithm  is  given  by  the  set  of  macro  statements: 

LOOP, NEWT,  10 

CTAN 

FORM 

PCG  .LINE.1. ,100. 

NEXT 

where  the  LOOP  for  the  NEWTon  step  is  to  be  executed  for  a  maximum  of  10  iterations  (the 
iteration  will  terminate  earlier  if  the  prescribed  tolerance  on  energy  is  met),  CTAN  indicates  a 
compressed  tangent  is  to  be  constructed  as  described  above,  and  PCG  indicates  that  the  Prccon 
ditioned  Conjugate  Gradient  algorithm  is  to  be  used  with  LINE  search  (omit  LINE  if  no  line 
search  is  desired  -  the  current  search  requires  several  evaluations  of  the  residual,  consequent h 
for  large  problems  may  he  time  consuming),  the  I.  indicates  that  the  value  ol  <■>  is  unity,  and 
100  is  the  maximum  number  of  PCG  iterations  to  be  used  to  solve  the  equations  Mictna 
lively,  the  commands. 
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LOOP, NEWT,  10 

n  an,-— ,1 

LANCM.INE.1.,100. 

NEXT 

may  be  used.  The  nonzero  value  on  the  CTAN  instruction  indicates  that  a  residual  >s  to  be 
computed  as  well  as  the  compressed  tangent  and  is  thus  equivalent  to  CTAN  and  FORM.  The 
use  of  LANC  instead  of  PCG  indicates  that  the  Lanczos  algorithm  is  to  be  employed. 


8.  Numerical  Kxamplcs 


The  preconditioned  conjugate  gradient  algorithm  (as  well  as  the  Lanczos  method)  has 
been  tested  on  a  series  of  test  problems.  In  order  to  assess  the  overall  efficiency  of  problem 
solution  we  report  both  the  computing  time  (for  a  VAX  11/780  computer  operating  under 
UNIX  4  2  BSD)  and  storage  requirements  for  the  coefficient  array.  In  our  test  problems  we 
include  two  and  three  dimensional  solids  subjected  to  both  static  and  dynamic  loading  states 

Kxample  I.  Two  dimensional  cantilever  structure. 

The  first  example  considered  is  a  cantilever  structure  with  two  holes  to  induce  added 
stress  gradient  ellects  The  model  consists  of  225  nodes  with  184  4-node  plane  elements,  see 
Figure  I  The  material  is  linear  elastic  and  utilizes  hl.MTOl  described  in  Chapter  24  ol  (13). 
For  this  problem  (as  well  as  all  subsequent  analyses)  we  perform  a  solution  using  the  direct 
solution  of  the  equations  as  well  as  the  PCCi  and  Lanczos  algorithms.  The  liming  and  perfor¬ 
mance  of  the  PCCi  and  Lanczos  methods  utilized  are  nearly  identical,  accordingly,  we  shall 
report  only  the  results  for  the  PCCi  algorithm.  The  essential  results  for  the  cantilever  struclure 
are  summarized  below. 

Model  Cantilever  type  structure 
(225  nodes  184  elements) 
profile  9990  Non-zero  terms  1162 

Static 

Direct 

total  time  16.77 

PCCi 

total  time  28  78  <  39  iterations  > 

Dynamic 
(5  time  Direct 
steps)  total  time  32.42 
PCCi 

total  time  77.52  •  24  to  27  iterations> 

These  results  are  much  as  expected  -  indicating  that  the  iterative  PCCi  algorithm  requires  more 
solution  ellort  (measured  in  CPU  time)  than  a  direct  solution.  The  only  redeeming  feature  for 
this  example  is  the  reduced  storage  requirements  for  the  stiffness  matrix  (3162  words  instead  of 


9940  words)  Accordingly,  if  one  had  access  to  a  very  small  computer  it  is  conceivable  that  the 
PCC»  algorithm  would  be  more  effective  since  it  could  greatly  reduce  the  number  of  calls  to 
backing  storage.  On  the  otherhand,  with  access  to  a  virtual  memory  machine  the  direct  solu¬ 


tion  is  to  be  preferred. 

Example  2.  Cylindrical  Structures 

As  a  second  example  we  consider  a  cylindrical  structure  subjected  to  end  loadings.  Two 
different  meshes  are  considered  to  illustrate  the  performance  of  the  PLM  algorithm  under  mesh 
refinements  The  material  is  again  linearly  elastic  and  both  static  and  dynamic  loadings  are  con¬ 
sidered.  The  first  mesh  consists  of  231  nodes  and  200  4-node  isoparametric  elements  (type 
PLMTOI).  while  the  refined  mesh  consists  of  4%  nodes  and  450  elements.  The  meshes  are 
shown  m  figures  2.  and  3.  Results  for  the  analyses  are  summarized  below 


Model:  Small  Cylinder  Structure 
(231  nodes  200  elements) 
profile  17485  Non-zero  terms:  3345 


Static 

Direct 

total  time:  22.50 

PC’Ci 

total  time  30.63  <  39  iterations> 


Dynamic 
(15  time  Direct 
steps)  total  time 


P(  (i 


83  32 


<.  25  to  27  iierations> 


total  time:  232.75 


Model:  Large  Cylinder  structure 
(496  nodes  450  elements) 
profile  57280  Non-zero  terms:  7570 

Static 

Direct 

total  time:  71.82 

PCG 

total  time:  149.48  <  140  iterations> 

Dynamic 

(5  lime  Direct 

steps)  total  time:  119.73 

PCG 

total  time:  230.37  <31  to  39  iterations> 

Once  again  the  direct  solution  is  more  efficient  in  CPU,  however  the  storage  requirements  for 
the  PLM  (or  PCG)  method  are  significantly  less  than  the  direct  method.  Note  that  the  number 
of  terms  is  almost  directly  proportional  to  the  number  of  nodes  (indeed  for  a  regular  mesh  of 
4-node  quadrilateral  elements  the  number  of  nonzero  terms  in  any  column  is  10  or  less), 
whereas  for  the  direct  method  the  number  of  terms  within  the  nonzero  profile  of  the  matrix  is 
almost  proportional  to  the  number  of  nodes  squared!  The  other  significant  fact  in  this  example 
is  the  number  of  iterations  required  to  solve  the  dynamic  problem  is  significantly  less  than  that 
required  for  the  static  loading.  Furthermore,  for  the  dynamic  loading  the  number  of  iterations 
required  to  solve  the  problem  increases  very  little  with  increased  problem  size. 

Example  3.  Three  Dimensional  Structures 

In  order  to  assess  the  performance  of  the  PCG  algorithm  on  three  dimensional  problems 
we  have  considered  the  loading  on  a  compact  block  of  8-node  brick  isoparametric  elements 
Two  different  meshes  with  linear  elastic  material  properties  have  been  considered.  The  first 
mesh  consists  of  64  elements  which  are  arranged  in  a  regular  cube  with  4  elements  on  a  side 
1  he  mesh  has  12s  nodes  with  21795  words  required  to  store  the  nonzero  profile  for  a  direct 
solution  and  only  7455  words  required  for  the  PCG  method.  For  the  dynamic  loading  case,  this 
problem  produces  the  first  PCG  results  which  are  more  efficient  than  a  direct  solution.  The 
results  are  summarized  below: 


Model  4X4X4  solid  structure 
<  125  nodes  64  elements) 
profile  21795  Non-zero  terms:  7455 


Static 

Direct 

total  time:  46.65 

PCG 

total  time:  201.49  <  187  iterations> 


Dynamic 
(4  time  Direct 
steps)  total  time:  93.18 

PCCi 

total  time:  79.13  <25  to  27  iterations> 

In  order  to  assess  the  improvement  in  performance  we  constructed  a  larger  problem  by  subdi¬ 
viding  the  mesh  to  form  a  cube  with  8  elements  on  each  edge.  Accordingly,  the  mesh  now 
contains  512  8-node  brick  elements  with  a  total  of  729  nodes.  The  nonzero  profile  increases 
dramatically  to  469,071  words  whereas,  as  before,  the  number  of  nonzero  terms  in  the 
compressed  profile  only  increases  proportionally  to  the  number  of  nodes  to  60,903  words.  The 
ratio  of  solution  times  for  the  dynamic  loading  case  increases  even  more  for  this  case,  as  sum¬ 
marized  below,  moreover,  even  the  static  loading  case  now  requires  less  CPU  for  the  PC’Ci 
method  than  that  of  the  direct  solution. 


Model:  8X8X8  solid  structure 
(729  nodes  512  elements) 
profile  469071  Non-zero  terms:  60903 


Static 

Direct 

total  time:  1585.23 

PCG 

total  lime: 1 145.29  <  130  iteralions> 


Dynamic 
(4  time  Direct 

steps)  total  time:2070.32 
PCG 

total  lime  320.83  <  5  and  6  iterations> 


This  example  illustrates  the  type  of  improvements  which  should  be  attained  for  all  large  three 


dimensional  applications  I  hc  size  of  problem  we  have  considered  is  guile  small  (indeed  even 
the  largest  mesh  we  could  consider  is  only  marginally  acceptable  lor  simulating  very  simple 
geometries)  and  is  limited  primarily  by  the  fact  that  we  utilized  a  VAX  11/780  computer.  In 
double  precision  arithmetic  we  required  over  4  megabytes  of  dimensioned  memory  to  solve  the 
problem.  We  fully  anticipate  that  applications  to  larger  problems  on  faster  ane  larger  computers 
can  achieve  the  same  level  of  improvement  we  have  indicated  here. 

Example  4.  Nonlinear  Material  Response  -  Two  Dimensional  Application. 

In  order  to  test  the  performance  of  the  PLM  algorithm  in  a  nonlinear  application,  we  con¬ 
sidered  the  elastic-plastic  static  response  of  a  plane  strain  strip  with  a  hole.  The  mesh  is  shown 
in  Figure  4  and  the  spread  of  plastic  zone  at  different  load  steps  in  Figure  5,  The  problem  was 
solved  using  both  direct  solution  and  the  PLM  method  and  utilized  the  consistent  tangent  for¬ 
mulation  developed  in  [121.  This  formulation  ensures  a  quadratic  asymptotic  rate  of  conver¬ 
gence  when  used  with  a  full  Newton  method.  The  overall  solution  time  for  the  PLM  method 
was  gt eater  than  the  direct  solution,  in  accordance  with  results  obtained  for  Examples  1  and  2. 
The  PLM  algorithm  performed  well,  however,  and  showed  no  loss  in  performance  with 
increased  plastic  deformations.  Accordingly,  we  fully  expect  that  the  solution  of  nonlinear 


three  dimensional  problems  will  be  more  efficient  with  the  PLM  method  than  a  solution 
achieved  using  a  direct  solution  of  the  algebraic  equations 
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9.  C  omparison  of  (ilohal  and  Element  by  Element  Preconditioned  Methods 

In  order  to  indicate  the  number  of  operations  in  an  element-by-element  (EXE)  method 
compared  to  the  global  preconditioned  form  with  compressed  storage  of  the  array  we  have 
constructed  a  table  to  indicate  the  number  of  operations  in  a  single  iteration  of  each  meth¬ 
od.  For  the  element-by-element  method  we  assume  a  second  order  accurate  double  pass 
method  (e.g.,  see  |8]).  The  results  are  summarized  in  the  table  for  Examples  1,  2,  and  3 
cited  above. 


I  Example 

! 

Mesh 

Operations  per  Iteration  | 

I 

' 

i 

Firm 

PCCi/PLM 

FXK 

i 

r  ” 

184 

12,648 

23,552 

i 

i  2 

200 

13,380 

25,600 

t 

!  2 

450 

30,280 

57,600 

|  3 

64 

29,820 

73,728 

1 

1  3 

1 

512 

243,612 

589.824 

The  difference  in  the  number  of  operations  is  due  to  the  fact  that  each  degree-of-freedom  m  a 
mesh  is  associated  with  more  than  one  element.  Indeed,  on  the  average,  the  above  table  indi¬ 
cates  that  there  is  a  savings  in  number  of  operations  by  a  ratio  of  about  1.8  to  2.5  tor  the 
PCCi/PLM  algorithms  in  comparison  with  an  element-by-element  algorithm  Thus,  an 
element-by-element  method  must  converge  in  about  half  as  many  steps  in  order  to  be  as 
efficient  as  the  PCCi/PLM  methods.  Our  previous  experience  indicated  that  element  precondi¬ 
tioning  never  converged  in  fewer  steps  that  the  global  preconditioning  method;  consequently, 
we  believe  that  the  current  implementation  offers  considerable  savings  over  element-by- 
elemcnt  methods.  The  final  proof  of  this  assertion  must,  however,  await  considerable  numeri¬ 
cal  testing  of  various  implementations  lor  iterative  methods. 
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Appendix  A.  Listing  of  C  ompact  Stiffness  Construction. 

The  lull  description  of  the  algorithm  to  construct  the  compact  storage  of  the  stiffness  is 
described  in:  "An  Algorithm  for  Assembly  of  Stiffness  Matrices  into  a  Compacted  Data  Struc¬ 
ture."  by  B  Nour-Omid  and  R.  L.  Taylor,  Report  No.  IJCB/SESM  84/06,  Structural  Engineer¬ 
ing  and  Structural  Mechanics,  University  of  California,  Berkeley,  May  1984  The  listing  fol¬ 
lows: 

SUBROUTINE  ELCNT ( NUMNP , NUMEL .NEN.NENI . IX , IC) 

DIMENSION  I X( NEN 1 . 1 ) , IC( t ) 

C 

C  INPUT 

C  NUMNP 

C  NUMEL 

C  NEN 

C  NEN  1 

C  IX 

C 

C  OUTPUT 

C  IC 

C 
C 
C 

C  COUNT  THE  NUMBER  OF  ELEMENTS  EACH  NODE  BELONGS  TO 

C 

CALL  I ZERO( IC  NUMNP ) 

DO  110  N  =  1 . NUMEL 
DO  100  J  =  1 , NEN 
I  .-  I  X  (  J  ,  N  ) 

I F  (  I  GT  0  )  IC(  I)  ■•=  IC(  I  )  +  1 
100  CONTINUE 

110  CONTI  NUE 

C 

C  SET  UP  POINTERS 

C 

DO  1  2  0  I  :=  2  .  NUMNP 

I C (  I  )  -  I C (  I  )  +  I  C(  I  -  1  ) 

120  CONTINUE 

C 

RETURN 
END 


SUBROUTINE  CASSEM(  D ,  A  ,  B  ,  S  .  P,  JCOLE  .  IRCAV,  LD.  ID,  NST.NEL  ,  AFL  .  BFL  ) 
IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

LOGICAL  AFL, BFL 

DIMENSION  D(l),A(l).B(l) , S(NST,  1  ) , P(  1  )  , JCOLE ( 1  )  ,  I  RCW(  I  )  , LD(  1  ) 
1  . ID( t ) 

C 

C  COMPACT  ASSEMBLY  OF  PROFILE  MATRIX 

r 

DO  200  J  I , NEL 
N  LD(J) 

IF  (  AFL  AND  N  GT  1  )  THEN 

DO  I  SO  I  1  ,  NEL 

K  l,D(  I  ) 

IF  (  K  GT  0  AND  K  LT  N  )  THEN 

INZ  ••=  I  NZ  A  (  JCOLE  (  N  -  1  )  +  1  .  JCOLE  (  N )  ,  I  RCAV,  K  ) 

A(  INZ  )  -=  A(  INZ  )  f  S (  I  ,  J  ) 

END  I  F 

ISO  CONTINUE 

END  IF 

C  ASSEMBLE  THE  DIAGONAL 

IF  (  N  GE  1  )  THEN 

IF  (  AFL  )  D(N)  -  D(N)  +  S ( J , J ) 

C  ASSEMBLE  THE  LOAD  IF  NECESSARY 

IF  (  BFL  )  B  (  N )  B  (  N )  ■+  P  (  J  ) 

END  IF 

200  CONTINUE 
RETURN 
END 


TOTAL  NO  OF  NODES  IN  THE  MESH 
TOTAL  NO  OF  ELEMENTS  IN  THE  MESH 
MAX  NO  OF  NODES  PER  ELEMENT 
DIMENSION  OF  IX  ARRAY 
ELEMENT  CONNECT  I  V  I TY  ARRAY 


ARRAY  OF  LENGTH  NUMNP  IT  FIRST  HOLDS  THE  ELEMENT  DEGREE 
OF  EACH  NODE.  THEN  BECOMES  A  POINTER  FOR  AN  ARRAY  THAT 
CONTAINS  THE  SET  OF  ELEMENTS  CONNECTED  TO  EACH  NODE 
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r. 


i 
•  .* 


c 

c 

c 

c 

c 

c . 

c 

c 

c 

c 

c 

c 


200 


:  1  0 


2  2  0 
230 
C 

c 

c 


r 

C 

C 


300 


3  1  0 
320 
330 


3  4  0 
360 


SUBROUTINE  COMPROf  NUMNP , NVMEL , NEN , NEN1  , NDF ,  IX,  ID, 1C.  IRCW,  IELC, 

1  JCOLE.KP) 

DIMENSION  I X(NENl , 1 ) , ID(NDF , 1 ) , IC( 1 ) , IRCW( 1 ) , IELC( I ) , JCOLE( I ) 

FOR  (NUMNP .NUMEL .NEN.NENl , IX, 1C)  SEE  SUBROUTINE  ELCNT 
INPUT 

NDF  NUMBER  OF  UNKNOWNS  AT  EACH  NODE 

ID  ACTIVE  UNKNOWNS  AT  EACH  NODE 

OUTPUT 

IELC  HOLDS  THE  SET  OF  EL EMENTS  CONNECTED  TO  EACH  NODE 
I  ROW  ROW  NUMBER  OF  EACH  NONZERO  IN  THE  STIFFNESS  MATRIX 
JCOLE  END  OF  ENTRIES  IN  I  ROW  FORM  A  GIVEN  COLUMN 

FIND  ELEMENTS  CONNECTED  TO  NODES 

CALL  I  ZERO  ( IELC,  I C( NUMNP ) ) 

DO  230  N  l  ,  NUMEL 
DO  2  2  0  J  =  I  , NEN 
I  =  I X( J , N) 

IF  (  I  GT  0  )  THEN 
KP  =  I C  (  I  ) 

IF  (  IELC(KP)  EQ  0  )  GO  TO  210 
KP  =  KP  -  1 
GO  TO  200 
IELC(KP)  -  N 
END  IF 
CONTINUE 
CONTINUE 

SET  UP  COMPRESSED  PROFILE  POINTERS 

KP  —  0 
NEP  =  I 

DO  360  I  =  1 . NUMNP 
NE  —  I C( I ) 

DO  3  4  0  II  -  I  ,  NDF 
NEQ  -  ID( I  I  ,  I ) 

IF  (  NEQ  GT  0  )  THEN 
JCOLE (NEQ)  -  KP 
KPO  =  KP  +  1 
IF  (  NEP  LE  NE  )  THEN 
DO  3  3  0  N  --  NEP ,  NE 
NN  -=  IELC(N) 

DO  3  20  J  --  I  ,  NEN 
K  I  X(  J  ,  NN ) 

DO  31 0  J  J  =  1 , NDF 
NEQ I  =  I D( J  J , K) 

IF  ( NEQJ  GE  NEQ  OR  NEQJ  LT  0 )  GO  TO  310 

CHECK  TO  SEE  IF  NODE  ALREADY  IN  LIST 

IF  (  KPO  LE  KP  )  THEN 
DO  300  KK  =  KPO , KP 

I F (  IRCW(KK)  EQ  NEQJ  )  GO  TO  310 
CONT I NUE 
END  IF 
KP  =  KP  +  1 
IRCW(KP)  ---  NEQJ 
CONTINUE 
CONTINUE 
CONTINUE 
JCOLE  (NEQ)  KP 
END  IF 
END  I  F 
CONTINUE 
NEP  —  NE  +  1 
CONTINUE 
RETURN 
END 


2, 

>.v 


INTEGER  FUNCTION  I NZA( Nl  , N2 , I  ROW, K ) 
DIMENSION  IRCW(l) 

FIND  THE  TERM  FOR  THE  ASSEMBLY 

DO  100  N  =  Nl  ,N2 

IF  (  IROW(N)  .  EQ  K  )  THEN 
INZA  =  N 
RETURN 
END  IF 
CONTINUE 

ERROR  IF  LOOP  EXITS 

STOP 

END 


SUBROUTINE  IZERO(IA.NN) 
DIMENSION  IA(NN) 

DO  100  N  =  1  ,  NN 
I A(N)  =  0 
CONTINUE 
RETURN 
END 


